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We study the in-plane and out-of-plane density ordering instabilities of quasi-two-dimensional fermionic po- 
lar molecules in single-layer and multi-layer configurations. We locate the soft modes by evaluating linear 
response functions within the conserving time-dependent Hartree-Fock (TDHF). The short-range exchange ef- 
fects are taken into account by solving the Bethe-Salpeter integral equation numerically. An instability phase 
diagram is calculated for both single-layer and multi-layer systems and the unstable wave-vector is indicated. 
In all cases, the in-plane density wave instability is found to precede the out-of-plane instability. The unstable 
wave-vector is found to be approximately twice the Fermi wave-vector of one of the subbands at a time and can 
change discontinuously as a function of density and dipolar interaction strength. In multi-layer configurations, 
we find a large enhancement of density wave instability driven by dilute quasiparticles in the first excited sub- 
band. Finally, we provide a simple qualitative description of the phase diagrams using a RPA-like approach. 
Compared to previous works done within the RPA approximation, we find that inclusion of exchange interac- 
tions stabilize the normal liquid phase further and increase the critical dipolar interaction strength corresponding 
to the onset of density- wave instability by over a factor of two. 



I. INTRODUCTION 



The field of ultracold atoms has witnessed a rapid progress 
in the past decade. Much of this experimental and theoreti- 
cal progress has been motivated by the prospect of realizing 
novel strongly correlated quantum phases and exploring the 
consequences of strong interactions CHll- One of the lat- 
est breakthroughs in this direction is the experimental real- 
ization of nearly quantum degenerate gases of fermionic polar 
molecules. By association of atoms via a Feshbach resonance 
to form deeply bound ultracold molecules a nearly de- 

generate gas of KRb polar molecules has been recently real- 
ized in their rotational and vibrational ground state 16 -101. 
The molecules can be polarized by the application of a dc 
electric field, resulting in strong dipole-dipole inter-molecular 
interactions. 

At the time this paper is being written, the coldest realized 
gas of polar fermionic molecules has a temperature of 1.4 
in the experiments of the group at JILA | 6 - 10], where Tp is 
the Fermi temperature. A major obstacle toward further evap- 
orative cooling of a large class of bi-alkali polar molecules 
is the existence of an energetically allowed two-body chem- 
ical reaction channel ifTTI . resulting in significant molecule 
losses in two-body scatterings. In a low temperature gas com- 
posed of a single hyperfine state, Fermi statistics blocks scat- 
terings in the 5'- wave channel and the majority of scatterings 
take place through the /7-wave channel. In a three-dimensional 
gas, the attractive head-to-tail dipolar interactions soften the 
/7-wave centrifugal barrier and increase the cross section of re- 
active collisions. The rate of chemical reactions can be effec- 
tively suppressed by loading the gas into a one-dimensional 
optical lattice (or trap) and aligning the dipoles perpendicu- 
lar to the formed quasi-two-dimensional layers, also known 
as pancakes. In such geometries, the incidence of head-to- 
tail scatterings is effectively suppressed due to the transverse 
confinement of the gas on one hand, and reinforcement of the 
/7-wave barrier due to repulsive side-by- side dipolar interac- 
tions on the other hand ||9l[T2l[T3l. Therefore, the preferred 
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FIG. 1. (Color online) (a) A multi-layer system of quasi-two- 
dimensional polar molecules. The dipoles are aligned perpendicu- 
lar to the x-y plane by the application of a strong dc electric field. 
The inter-layer spacing is d. The x-axis is perpendicular to the plane 
of the plot, (b) each quasi-two-dimensional layer is composed of 
multiple subbands, corresponding to the transverse wavefunctions of 
the trap. For symmetric lattice potentials, the subbands have a well- 
defined parity with respect to the reflection about the x-y plane. The 
trap confinement width, a±, is shown in the figure. 



geometry to study reactive polar molecules is in tightly con- 
fined two-dimensional layers. 

In such geometries, the energy levels of particles is 
quantized due to the transverse confining potential and each 
quasi-two-dimensional layer can be thought of as a collection 
of two-dimensional energy subbands (see Fig. [T]). Since 
higher subbands have a larger transverse spreading, it is 
expected that occupation of higher subbands will increase 
the rate of head-to-tail collisions and consequently, the 
molecule loss rate. However, it has been recently shown that 
the two-body chemical reactions will still be significantly 



suppressed even if the first few subbands are filled, due to 
Pauli blocking |13|. The occupation of higher subbands does 
not impose any difficulty on experiments with non-reactive 
species such as NaK, NaRb, NaCs, KCs, and RbCs |11J. 
The possibility of going beyond the single- subband limit 
opens a new window toward experimental and theoretical 
exploration of many-body physics of quasi-two-dimensional 
fermionic systems with anisotropic interactions. 

In contrast to the isotropic short-range interactions in ul- 
tracold atomic gases realized using a ^--wave Feshbach reso- 
nance the long-range and anisotropic nature of elec- 
tric dipole-dipole interactions in ultracold gases of polar 
molecules allows the experimental realization of a wider range 
of physical phenomena. In particular, the repulsive side- 
by-side dipole-dipole interactions in layered stacks of polar 
molecules can lead to spontaneous translational symmetry 
breaking and formation of density ordered phases for strong 
interactions. We define the ratio of the typical interaction 
over the kinetic energy, r^^, as a dimensionless measure of the 
strength of dipolar interaction: 



where D is the electric dipole moment of a single molecule, 
m is the molecular mass and n is the two-dimensional 
density. It is noticed that in contrast to the electron gas, the 
interaction energy is dominants at higher densities for fixed 
dipole moments and as a result, the density ordered phases 
are expected to appear at higher densities. 

Recently, Yamaguchi et ah 1 14| and Sun et al. fTS] have 
independently studied the density-wave (DW) instability in 
a strictly two-dimensional layer of polar molecules in the 
random phase approximation (RPA). At zero temperature and 
zero tilt angle of dipoles with respect to the confining 2D 
plane, their results indicate that the DW instability occurs for 
r^^ ^ 0.17. The former study treats the self-energy corrections 
within an approximate variational method and first-order 
perturbation theory |[T4ll . The second study neglects the self- 
energy corrections altogether, however, present a rigorous 
proof for the necessity of DW instability for strong enough 
interactions and predict the nature of the density ordered 
phases at different tilt angles |15|. Both studies neglect the 
exchange interaction effects beyond the cancellation of the 
5-- wave component of dipolar interaction in their calculations. 
Since the interactions need to be appreciably strong for the 
density ordered instabilities to occur and that the ordering 
wave- vector is in order of the Fermi wave- vector, we expect 
that inclusion of short-range exchange effects will result in a 
significant quantitative correction to the results of the cited 
works. 

The simplest self-consistent and conserving many-body 
approximation that respects the Fermi statistics is the time- 
dependent Hartree-Fock approximation (TDHF) |16|, also 
known as the generalized random phase approximation 
(GRPA) [TTtI. We have recently studied the band renor- 
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FIG. 2. (Color online) A schematic representation of the homoge- 
neous liquid and density ordered phases of a multi-layer system of 
quasi-two-dimensional polar molecules, (a) the liquid phase is char- 
acterized by uniform in-plane and out-of-plane density in each layer, 
(b) the ripplon phase is characterized by out-of-plane density mod- 
ulations and uniform in-plane projected density. The Z2 reflection 
symmetry may be broken (as shown here) if the mixing occurs be- 
tween subbands of even and odd parity. In that case, the energeti- 
cally favorable configuration corresponds to a tt phase shift between 
even and odd layers due to the inter-layer attraction, (c) the density- 
wave phase is characterized by broken in-plane translational sym- 
metry. Wigner crystals, striped and bubble phases are examples of 
density- wave phases, (d) the zigzag crystal phase is characterized by 
broken in-plane translation symmetry and presence of out-of-plane 
density modulations. 



malization and collective modes of quasi-two-dimensional 
polar molecules within the TDHF approximation flSl . In 
this paper, we study the density ordering instabilities of the 
liquid phase of quasi-two-dimensional polar molecules in 
single-layer and multi-layer configurations. Throughout this 
study, we assume that the layers are well- separated such that 
the inter-layer tunneling can be neglected. 

We consider instability toward two types of density ordered 
phases: in-plane density- wave phase and the ripplon phase 
(see Fig. [2]). The in-plane density-wave phase is characterized 
by broken in-plane translational symmetry (see Fig. [2]:). 
The ripplon phase is a reminiscent of the spin density-wave 
(SDW) phase of electronic systems 1 19], where quasiparticles 
of different subbands play the role of different spin states. 
Since quasiparticles in difference subbands have different 
transverse wavef unctions, their mixing results in out-of-plane 
density modulations (see Fig.[2j)). In the ripplon phase, the Z2 
reflection symmetry about the confining plane may be broken 
if the mixing occurs between subbands of even and odd parity. 
When the translational symmetry of the in-plane projected 
density is broken in addition to the presence of subband 
mixing, we denote the phase by zigzag phase (see Fig. |2]l). 
Quantum zigzag transition in one-dimensional ion chains 
has been a subject of active theoretical and experimental 
investigations in the past few years ll20K24ll . 

We adhere to the mode-softening paradigm of phase tran- 
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sition and look for the instabilities by monitoring various 
density-density response functions in the normal phase as the 
interaction and trap strengths are varied. The softening of a 
density-wave mode results in development of a sharp peak 
in its corresponding response function, which eventually pro- 
motes to a singularity when the liquid phase becomes unstable 
against density fluctuations. 

We remark that the TDHF approximation adopted in this 
study is essentially an analysis of Gaussian fluctuations about 
the mean-field liquid phase. In many cases, analyses at the 
mean-field level underestimate the stability of symmetric 
phases and predict transition to the symmetry broken phases 
too early. Inclusion of correlation effects usually enhance 
the stability of the unbroken phase beyond mean-field 
predictions. Thus, the instability criterion resulting from 
TDHF response functions is most likely only a signal for the 
formation of short-range correlations in the liquid phase, with 
the phase transition to the density ordered phase occurring at 
stronger interactions. Nevertheless, the mean-field stability 
analysis is a first step and an indispensable guide in construct- 
ing more elaborate approximations. We discuss this issue 
with more detail in Sec.lVll 

This paper is organized as follows: we review the micro- 
scopic model for single-layer and multi-layer configurations 
in Secim The TDHF formalism for multi-subband and multi- 
layer systems is described in S ec.[lll| and the instability phase 
diagrams are presented in Sec. |IV| An approximate RPA-like 
theory is developed in Sec.[V|using which the qualitative fea- 
tures of the obtained instability phase diagrams are explained. 



Finally, we discuss the results in Sec. VI and compare our 
results with the previous works. The derivation of analyti- 
cal expressions for the effective inter-layer and inter-subband 
dipolar interactions and the details of our numerical method 
are provided in the appendices. 



II. THE MICROSCOPIC MODEL 

We start our analysis by reviewing the microscopic model 
describing fermionic polar molecules of mass m, prepared in a 
single hyperfine state, and loaded in a one-dimensional optical 
lattice. For concreteness, we assume that the optical lattice is 
along the z-axis. Also, we assume throughout this paper that 
all of the dipoles are aligned perpendicular to the trap plane 
by the application of a strong external d.c. electric field (see 
Fig.[T]). We work within the units H = 2m = 1 unless these 
quantities appear explicitly. We also denote 3D and in-plane 
2D coordinates by r and x respectively. 

The second-quantized Hamiltonian describing the system 
can be written as: 



n 



(2) 

where ?/^(r) is the fermion annihilation operator and Mat.(^) 



and Vdip.(i* ~ 1*0 denote the optical potential and dipolar in- 
teraction respectively: 



^dip.(r) 



(|rp-3.^) 



(3) 
(4) 



where A is the wavelength of the optical lattice lasers. The gas 
in the optical lattice can be thought of as a stack of quasi-two- 
dimensional layers separated by a distance d = A/2. In prac- 
tice, there are a finite number of layers present in the sample. 
We denote the number of layers by Ni and assume a periodic 
boundary condition along z direction in order to eliminate the 
surface effects. We also denote the transverse size of the stack 
by L = Nid. 

The fermion operator can be conveniently expanded in 
Wannier's basis in z direction and plane- wave basis in x-y 
plane: 



oo Ni 



k a=l n=l 



(5) 



where A is the area of system in x — y planes, Wan{z) de- 
notes the Wannier's wavefunction of the band a, with its cen- 
ter shifted to n'th well of the optical lattice and Cna,k anni- 
hilates a fermion in layer n, subband a and with momentum 
k. We omit the limits in the summations over layer and sub- 
band indices in the rest of the paper for brevity. Plugging the 
expansion of the fermion operator into Eq. (|2]), we get: 

k mn a(3 ki,k2,q 



Q!/3;7A mn;rs 



(q)4 



mo;, ki+q r7,k2 — q sA,k2 n/3,ki ' 



(6) 



where: 



J™/ = J dzw^Jz) 



W0n{z), (7) 



and: 



X w^^{z') Wsx{z') Wnisiz) Vdip(r - r'). 

(8) 

In order to simplify the analysis, we assume that the opti- 
cal potential is deep and that its minima are well separated, so 
that we can neglect inter-layer tunneling effects. We call this 
limit the independent layers (IL) limit. The subsequent results 
presented in this paper are all within this limit. In the IL limit, 
the optical potential (Eq. [Sj can be expanded to quadratic or- 
der about its minima: 
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(9) 
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where entrap = (27r/A) y^2Fo/^ in the effective trap fre- 
quency of each well. The Wannier's wavef unctions can also 
be approximated by shifted harmonic oscillator wavefunc- 
tions: 



■■Ha{z/a^), 



(10) 



where Ha is the Hermite polynomial of degree a, and a± is 
the transverse spreading of the lowest subband which is re- 
lated to the parameters of the optical lattice as: 



hX 



(11) 



The conditions of the IL limit is met if the transverse spread- 
ing of the Wannier's wavefunctions is smaller than the inter- 
layer separation, i.e. a± <C d. For sinusoidal optical poten- 
tials, we get the explicit condition 4/i^/Av^2mVo 1- 

In the IL limit, the overlap between the wavefunctions of 
different layers is negligible and one can assume: 

^am{z) Wj3n{z) (X 5mni fo^ ^11 z G M, and m, n G Z. 

(12) 

It is straightforward to show using Eqs. (|7]), ^ and ( 12 ) that 
the one-body and two-body matrix elements appearing in the 
second-quantized Hamiltonian assume a much simpler form 
in this limit: 



J, 



v: 



Q!/3;7A 



al3 

(q) 



-'mn'-'rs 



(13) 
(14) 



where is the zero-point energy of a'th subband and 
is explicitly given by hujtrapio^ + 1/2) in the harmonic 
trap approximation described above. We note that a more 
explicit condition for the IL limit in the negligibility of the 
off-diagonal matrix elements of and in the 

layer indices. Intuitively, Eq. ([13]) and ([14]) imply the absence 
of inter-layer tunneling and inter-layer exchange interactions, 
respectively. It is easy to see that the layer index remains 
a good quantum number in the IL limit in the presence of 
interactions and this salient feature simplifies the study of 
multi-layer systems to a great degree. 

In the next two subsections, we explore the effective micro- 
scopic models for single-layer systems (Ni = 1) and multi- 
layer cases (Ni > 1) in some more detail and briefly discuss 
the features of the normal liquid phase in each case within the 
Hartree-Fock approximation. 



A. Single-layer systems 

The absence of inter-layer attractive interactions in a single- 
layer system makes it an ideal starting point for the study 
of the more complicated case of a multi-layer configuration. 
The physics of single-layer systems is essentially governed by 
intra-layer repulsive interactions. From an experimental point 



of view, this limit is achieved either by selectively removing 
particles from an optical lattice in order to get a single pan- 
cake, or by utilizing a strong optical trap instead of an optical 
lattice. The Hamiltonian takes the following form in this limit: 



Q!,k ki ,k2 ,q Q!/3;7A 

^a(3;-fX{^) c]^,ki+q ^7,k2-q^A,k2 ^/3,ki ' 



(15) 

where Va(3--fx{(l) = "^a^-^xi^) intra-layer interaction 

matrix elements. A generating function and explicit expres- 
sions for Vaf3;-fx{^) can be found in our earlier paper |18|. 
For trap potentials which are symmetric about their center, 
the effective inter-subband interaction matrix elements 
conserve the net subband parity of the scattering particles, i.e. 
Vc/3;7A(q) /0ifa + /3 + 7 + A = 0(mod2) (13. 

In an earlier paper, we have studied the self-energy correc- 
tions of the single-layer system in the normal liquid phase in 
Ref. 1 18] in the self-consistent Hartree-Fock approximation. 
We do not repeat the analysis here and just mention that be- 
sides the usual Hatree-Fock band renormalization, one also 
finds that the non-interacting subband indices do not remain 
good quantum numbers in the presence of interactions. A 
well-defined subband index can still be found after applying 
an orthogonal transformation that diagonalizes the Hartree- 
Fock decoupled Hamiltonian. More explicitly, one can define 
a set of Hartree-Fock fermion operators as: 



such that: 



k,CK ^k,a' 



(16) 



(17) 



Q!,k 



where eQ,(k) are the renormalized energy dispersions of the 
Hartree-Fock subbands and 1-L§[ is the Hartree-Fock decou- 
pled Hamiltonian of a single-layer system. The orthogonal 



transformations appearing in Eq. (16), /7^Q,(k), as well as 



the renormalized dispersions, ec^(k), are found by solving the 
Hatree-Fock equations (see Ref. 1 18 1 for details). The temper- 
ature Green's function for Hatree-Fock quasiparticles can be 



read directly from Eq. ( [17] ): 



Jo 



(18) 



where pg^ 



-fj'J^) /ZsL is the grand canonical equi- 



librium density matrix and ^/i(k) = e^(k) — /i. The Green's 
function in the original non-interacting basis can also be found 



using the inverse of the transformation given in Eqs. ( 16 ): 

/7^A(k)/7,A(k) 



■6(k) 



(19) 



The above expression for the Green's function is found to be 
useful in evaluating frequency summations later. 
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B. Multi-layer systems 

The physics of multi-layer systems is governed by both 
intra-layer and inter-layer interactions. As we will see later, 
the interplay of these forces will modify the density- wave in- 
stability of the system to a great degree. The Hamiltonian 
takes the following form in this limit: 



k a,m ki,k2,q 

/ ^ / ^ ^a(3;-fX ^rna,ki+q ^r7,k2-q 

Q!/3;7A mr 

XCrA,k2^m/3,kr (20) 

In contrast to the interaction of particles within each layer, 
the inter- subband interaction of particles across the layers 
violate the parity conservation due to the absence of re- 
flection symmetry, i.e. V^^~^^(q) can still be non-zero if 
a + /3 + 7 + A = l (mod 2), provided that m ^ r. Explicit 
expressions for V^^.'^]^(q) are provided in Appendix [a[ 

As mentioned in Sec. |ll| one simplifying aspect of the IL 
limit is the conservation of the layer indices in the scattering 



processes (see Eq.[T2]). As a consequence, the normal liquid 
solution of multi-layer systems in the Hartree-Fock approx- 
imation closely resembles that of single-layer systems. The 
Green's function is found by solving the following Dyson's 
equation: 



a 


a 


a 






m,q I 














/3 


/3 


f3 



m' ,k' 
vwwvl 1 



m.k' 



(21) 

where m and m' are layer indices, q is the momentum transfer, 
the greek letters denote subband indices and thin and thick 
lines denote bare and dressed Green's functions. The above 
diagrammatic equation yields: 



X ^M^;m(q)^^/3;m(q, ^^n), 



(22) 



where the non-interacting Green's function, ^^^.^(q, icj^), 
is: 



(23) 



and the proper self-energy matrix i;^^.^(q) is defined as: 

-VU'l-yM - k')^w] Qxr,m'{^', ico'J. (24) 



In the homogeneous normal liquid phase, the layers are iden- 
tical and the Green's functions and self-energy matrices are 
independent of the layer indices. As a result, the Hatree-Fock 
equation for multi-layer systems in the IL limit has an identi- 
cal structure to that of single-layer systems, however, with ad- 
ditional contributions coming from direct inter-layer interac- 
tions. Thus, the numerical method described in Ref. flE] can 
be identically applied to obtain the renormalized subbands of 
multi-layer systems as well. We refer the reader to Ref. iTSl 
for computational details and suffice it to mention that like 
single-layer systems, one can again find an orthogonal trans- 
formation of the bare fermion operators that diagonalizes the 
Hartree-Fock-decoupled Hamiltonian. More explicitly, one 
can define Hartree-Fock quasiparticle operators as: 



(25) 



such that: 



<l= E e~"(k)4a,k5™«,k, (26) 

m,Q!,k 

where (k) is the renormalized energy dispersion of subband 
a, UafdO^) is an orthogonal transformation and the 
Hartree-Fock decoupled Hamiltonian of a multi-layer system. 
Note that and [/^^ are the same for all layers. The temper- 
ature Green's function can be directly read from Eq. ([26]) and 
is formally identical to Eq. ([18]) and ( [T9| ) respectively. The ef- 
fects of direct inter-layer interactions are implicitly included 
in the renormalized subband dispersions and orthogonal trans- 
formations. 



III. EVALUATING THE RESPONSE FUNCTIONS IN THE 
TDHF APPROXIMATION AND LOCATING THE 
INSTABILITIES OF THE LIQUID PHASE 

The instabilities of the liquid phase can be located by cal- 
culating various static response functions in the liquid phase. 
We investigate the instability of the liquid phase toward in- 
plane and out-of-plane (ripplon) density- wave orders. As a 
first step, we define the order parameters and their correspond- 
ing response functions in more detail in the next subsections. 



A. Order parameters 

We define the in-plane projected density operator of layer 
m at in-plane coordinate x as: 

n{m+l/2)d 
J{m-l/2)d 

p{m-\-l/2)d 

= / dzw;^^{z)Wm^a^{z) 



-i(k-k')-x t 



EE^ "'~ci,„,,,a 



(k-k')-x^t 

mQ!,k'"mQ;,k' ' 



(27) 



Oi k,k' 
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where we have adopted the IL approximation in the last Hne. 
In the normal phase, (pm(x)) is constant and independent 
of X. The in-plane density- wave instability is characterized 
by appearance of (quasi-)periodic spatial modulations in 
{pm{^)) (see Fig.|2|:). 

We define the aP-ripplon operator of layer m at in-plane 
coordinate x as: 

^ J{ra-l/2)d ^ ^ 

k,k' 

(28) 

Again, we have adopted the IL approximation in the last 
line. In the normal phase, {S^{^)) = for a ^ /3. The 
a/3-ripplon instability is characterized by growth of (quasi- 
)periodic spatial modulations in {S^^{x.)) and absence of 
any instability in the in-plane projected density (see Fig.|2])). 
When both in-plane and out-of-plane symmetries are broken, 
we refer to the case as the zigzag instability (see Fig.[2]l). 

B. Evaluation of the response functions 

We evaluate the response functions in the imaginary time 
formalism and find the real-time response functions by an- 
alytic continuation. The imaginary-time in-plane projected 
density-density response function is defined as: 

X^r"'^(xr;xV) = -Tr{pGT. [p^(xr)p^, (xV)]} , 

(29) 

where pc = e~^^^~^-^^ /Zq is the grand canonical weight- 
ing operator. The imagiary-time ripplon-ripplon response 
function is defined as: 

= -Tr{pGT.K^(xr)^::,e(xV)]}. 

(30) 

The space and time translation invariance of the Hamiltonian 
and the normal phase imply that the response functions are 
functions of x — x' and r — r' and therefore, it is easier to 
express them in terms of transferred momentum q and bosonic 
Matsubara frequencies ivn- Also, both of the density-desity 
and ripplon-ripplon response functions can be expressed in 
terms of polarization insertions as follows: 



and: 



(m—m')/ I I 

Xa/3 '(xt;xt 



a,l3 



(31) 



and: 



(m — m')/ ■ \ 1 



n 



(m—m') 

a(3;a(3 



/ • \ I -|-r(m— m') / . \ 



ni™;;?^ ^ (q. ^'^n) + njj™.^^ ^ (q, w^) 



(32) 



where the polarization insertion is defined as: 

nS-r'^(q--») = ^EnS;;r'^(q'--k,k'), (33) 



^0 

(cL,k+q(^)c„.;3,k(^)4'^,k'-q(0)c™,A,k'(0)}^^^; (34) 
Only diagrams with connected external vertices must be 



considered in Eq. ( 34 ). 



The TDHF approximation for the polarization insertion 
amounts to summing ladder and ring diagrams to all or- 
ders CSIITtI. Although we are only interested in the static 
limit in this study, i.e. ^0+, we will only take this limit 

at the end of the derivation for generality. A typical term con- 
tributing to n[^^~^ ^ in the TDHF approximation consists of 
one or more bubble diagrams, possibly with ladder-type ver- 
tex corrections, connected to each other by interaction lines: 

j^(m-m') 



Q!/3;7A 




m'X 



m',7 



+ ... (35) 

Since the layer index in conserved on each interaction vertex, 
it is easy to see that the particle and hole lines appearing in 
an irreducible polarization diagram (bubble) carry the same 
layer index. Thus, the vertex corrections are only due to the 
intra-layer interactions. The homogeneity of the normal phase 
also implies that the bubble diagrams are independent of the 
layer indices. Thus, we can carry out the summation in two 
steps: first, we evaluate the irreducible polarizations by sum- 
ming the ladder-type vertex corrections to all orders. Next, we 
calculate the full polarization by connecting the bubbles with 
interaction lines. 



Let n 



Q!/3;7A 



be the irreducible intera-layer particle-hole 
propagator with ladder-like interactions summed to all orders. 
^a(3;-fX found by solving the following Bethe-Salpeter 

equation: 



a/3;'y A 



/5 v k'-\-q,a k2,X 



ki-\-q 



^ki,k2 



ki 



ki-\-q 
ki 



pcr;-yX 



(36) 



k,k' 



a 'J a fi k' ,p k^—qn 

The diagrammatic equation yields the following integral equa- 
tion: 

n«/3;7A(q.^^n;ki) = Ii^^l.^^{\^i,q^,iVn) 

X n*,.^;,(q,iz/„;k'), (37) 
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where we have summed both sides over k2 . Summation over 
repeated indices is assumed throughout. n[^^j.^^(q, i^n'i k) is 
the bare particle-hole propagator: 



Pi 



Uaa' (k) U^^> (k + q) U^^> (k) Uxx' (k + q) 



(^'k+q,/3' - ^k,a' 



(38) 



The irreducible polarization diagram, n*^.^;^(q, ii^^), is 
found by summing n^^.^;^(q, ii^^; ki, over ki and 
k2. The summation over k2 is trivial and is already 



done in Eq. (37). The summation over ki, however, may 
only be done once the solution of the integral equation is 
known. We solve the integral equation numerically. The 
details of the numerical procedure is provided in Appendix. [B] 

Once n^^.^;^(q, iun) is evaluated, the full polarization can 
be easily obtained by summing the ring diagrams to all orders. 
We note that the interaction lines connecting the irreducible 
polarizations may have vertices belonging to different layers 



(see Eq. 35 ). The following Dyson's equation yields the sum- 
mation ring diagrams to all orders: 



n 



{m—m') 
Q!/3;7A 



(q, i^i 




where the blank and filled circles denote irreducible and full 
polarizations, respectively. The above diagrammatic equation 
yields the following linear system of equations: 

^^a(3;-fX ~ ^rnm'i-^afdi-fX 
Ni-1 

n= — Ni-\-l ^iv\p(j 

where we have dropped the common argument {q_^iun) for 
brevity. Since we assumed periodic boundary conditions 

along the z-axis, n[^^~^ ^ is a periodic function of m — m' 
and the above equation can be most conveniently solved by 
going from the layer index to transverse momentum represen- 
tation. We define: 



Ka'xi'l^i'^n), (41) 



where kz = 27rn/L for n = 0, 1, . . . , A^^ — 1. Plugging 
Eq. ^ into Eq. (gO]), we get: 



n 



Kp;-yX + E n^ftM. '^i'^a K'n^ ' (^2) 



where: 



vjnfU^)= E (43) 

n=-Ni + l 

The transverse modes with different kz are decoupled in 
Eq. ( [42| ) and the problem reduces to solving a linear system 
in the subband indices for each kz. The response functions 
can also be expressed conveniently in the transverse momen- 



tum basis using Eqs. (31) and ( 32 ): 



Xdd 



(44) 



a, (3 



nSiJq.^^n) + n^;Y/5jq,z^0 



(45) 



Before embarking on evaluating the response function us- 
ing the described formalism, we find it worthwhile to briefly 
study the direct consequences of the coupling between in- 
plane and out-of-plane modes. Understanding the coupling 
between various density ordering modes guides us in pre- 
dicting which modes go unstable simultaneously and which 
modes may remain stable once the liquid phase becomes un- 
stable. 

It is straightforward to establish that all in-plane density 
fluctuations (corresponding to polarization diagrams such as 
noo;oo, noo;ii, IIii;!!, etc.) arc coupled to each other due to 
the existence of interaction matrix elements Voo;ii and such. 
Therefore, the in-plane density wave modes go unstable to- 
gether and contribute to the formation of an inhomogeneous 
case. In particular, coexistence of liquid phase in one subband 
and a density ordered phase in another subband is impossible. 



As mentioned in Sec. |II A| the inter- subband interactions 
conserve the net of parity of the interacting quasiparticles in 
the single-layer case. As a consequence, there is no coupling 
between in-plane density fluctuations and odd ripplons (cor- 
responding to polarization diagrams such as IIoijoi, Hoijio, 
etc) due to the absence of interaction matrix elements Voo;Oi 
and such. For instance, starting from the Ni phase, it is pos- 
sible to reach a density ordered phase with no accompanying 
out-of-plane order. 

In multi-layer systems (Ni > 1), the situation can be dif- 
ferent. As mentioned in Sec. |IIB| the inter-subband interac- 
tions between quasiparticles of different layers violate the par- 
ity conservation. Using the results of Appendix [A] one easily 
finds that the parity violating interaction matrix elements are 



V 



{m—m') 



q;/3;7A 



(q) 



odd under the inversion of layer indices, i.e. 

(-l)"" ^Lt^r^ (q)' where P = (a + /3 + 7 + A) mod 2. Us- 
ing this property, the inter-layer interactions in the transverse 
momentum basis, Eq. ([43]), can be expressed in a more useful 
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form: 

Ni-l 

V„ft^A(q) + 2 ^ cos(fc,nd)vi;)^^;,(q) 

n=l 

y(fc.) ( ifP = 0, 

-2i ^ sin(fc,nd)vl";;^;,(q) 

n=l 

if P = 1, 

(46) 

Clearly, the parity violating matrix elements (P = 1) are 
non- vanishing only if kz 7^ 0. Therefore, in multi-layer 
systems, density waves and odd ripplons are coupled at finite 
transverse momenta. Thus, if the first mode that goes unstable 
has a finite transverse momentum, the resulting ordered phase 
breaks both in-plane translation symmetry and Z2 reflection 
symmetry. 

We note that once the leading instability is found, the study 
of subsequent instabilities must be done with a word of cau- 
tion. The leading instability modifies the initial state, either 
by producing short range correlations or breaking a symmetry. 
Even if the new state can be described well at the mean-field 
level, the Green's function and the response functions must be 
recalculated in the new state. This requirement in turn modi- 
fies the nature and/or order of the subsequent instabilities. We 
only study the leading instability of the liquid phase in this 
paper and leave the study of subsequent transitions within the 
density-ordered phase for future works. 



IV. THE MEAN-FIELD INSTABILITY DIAGRAM OF THE 
LIQUID PHASE 

Due to the complexity of the formalism described in 
the previous section, obtaining analytical expressions for 
the response functions in the TDHF approximation is a 
formidable task without resorting to further approximations. 
The most involved part of the calculation is solving the 
Bethe-Salpeter integral equation that represents the effects 
of intra-layer exchange interactions. Here, we present the 
results obtained by exact numerical calculations based on the 
procedure outlined in the previous section. The numerical 
procedure is described in Appendix [B] in detail. Later, we will 
employ further simplifying approximations in order to obtain 
tractable analytical expressions that guide us in interpreting 
the results. 

We are interested in finding the boundary of the stability of 
the liquid phase and the characteristics of the mode that drives 
the instability, as a function of tunable parameters of the sys- 
tem. For a fixed number of layers Ni , inter-layer separation d 
and temperature T, there remains two tunable dimensionless 
parameters: the dipolar interaction strength, rd (see Eq. [T]), 
and the ratio of the transverse confinement width and the 
mean inter-particle distance, y/na^. The limits ^/na±_ <C 1, 
■\/naj_ ~ 1 and ^/naj_ 1 correspond to the two-, quasi- 
two- and three-dimensional regimes respectively. The IL limit 



is achieved for a±/d <C 1. In the study of multi-layer sys- 
tems, we restrict our parameters to ^/nd > 1, i.e. the high 
density limit (with respective to the layer spacing), in which 
both IL limit and quasi-two-dimensionality can be approxi- 
mately achieved. 

We locate the instability boundaries of the liquid phase us- 
ing a divide and conquer method. For each a_\_ , we first locate 
Td^L and Td^u such that all response functions are regular and 
smooth for rd^h and at least one mode is unstable at Vd^u. The 
instability appears as a zero crossing in the inverse of some 
response function. Once a lower and upper limit is found for 
the critical rd, the exact location of the phase boundary is de- 
termined by successive bisection of this interval. 

In order to simplify our analysis, we confine our atten- 
tion to the low temperatures, where thermal fluctuations are 
negligible compared to the quantum fluctuations. We set 
T = 0.02 T^, where T^^ = 2Tinh^/mkB is the Fermi tem- 
perature of a two-dimensional free Fermi gas at the same den- 
sity. We will later show that the chosen small temperature 
is high enough to suppress the inter-layer superfluid transi- 
tion 12511271 in all of the studied multi-layer configurations. 



A. Instabilities of single-layer systems 

We have studied the properties of the liquid phase of 
single-layer systems in an earlier paper 1 18|. In brief, when 
>/na±_ <C 1, the energy gap between the subbands is much 
larger than the Fermi energy and the system is effectively 
two-dimensional, i.e. only the lowest subband (a = 0) is 
filled. Upon relaxing the trap, i.e. y/na\_ ^ 1, the subband 
gap is reduced and higher subbands will be filled. We denote 
a normal liquid phase having up to j'th subband filled by Nj. 
The Fermi surface of a system in the Nj phase consists of j + 1 
circles, characterized by their radii kp^o, kp,!, ... , kpj- 
In analogy to quasi-two-dimensional electron gas, we expect 
to get j + 1 peaks in static density-density response function 
vs. momentum q at q ^ 2/ci?^07 Q ^ "^kp^i^ ... , q ^ ^/c^^, 
corresponding to softened particle-hole excitations arising 
from opposite poles of the Fermi surface of each subband. 
We also expect to get a single peak at q ^ kp^a ^ kp^is in 
the a/3-ripplon response function, again, analogous to SDW 
softening in electron gas |[T9l . 

Fig. |3] shows the static density-density response function in 
No (top plot) and Ni (middle and bottom plots) phases. In 
the No phase, we only get a single peak corresponding to the 
softened density- wave mode ai q ^ ^kp^o- The middle and 
bottom plots (Ni) correspond to low and high population of 
the first excited subband. It is noticed that in the middle plot, 
\hQ q ^ 2kp^Q mode is more enhanced compared io q ^ 2kp^i 
mode. The scenario is reversed, however, as the population 
of the first subband is increased beyond a certain threshold. 
Thus, we generally expect q ^ 2kp^o to be the first mode to 
go unstable in the No phase, while we expect a switching from 
q ^ 2kp^o to q ^ 2kp^i in the Ni phase. 

Fig. [4] shows the static 01-ripplon response function for the 
same configurations as in Fig.|3] A slight enhancement of the 



9 







I. . 




- . . , 1 . . . 


1 , , , 





n-V2|q| 



FIG. 3. Static density-density response function of a single-layer 
system in the normal phase. Xdd = 27r/i^Xdd/^ and the dashed 
lines denote q = 2k fj, j = 0, 1. (a) No phase [rd = 1.0, yjna^ 
= 0.15] (b) Ni phase [r^ = 1.35, ^/na± = 0.25] (c) Ni phase [r^ = 
1.35, ^/na±=035l 
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FIG. 4. Static 01-ripplon response function of a single-layer system 
in the normal phase. Xoi = ^nh^xoi / m and the dashed lines denote 
q = kF,o-\-kF,i- (a) No phase [r^ = 1.0, y/na± =0.15] (b) Ni phase 
[rd = 1.35, = 0.25] (c) Ni phase [rd = 1.35, V^a± = 0.35]. 



01-ripplon mode at = A:f,o + ^f,i is noticed in Ni phase, 
however, the peaks are less pronounced than the peaks of the 
density-density response function. This result can be under- 
stood in light of the stronger intra- subband vs. inter- subband 
repulsion, the latter being weaker due to contributions from at- 
tractive head-to-tail dipole-dipole interactions. Therefore, we 
generally expect the density-wave instability to preceed the 
ripplon instability. 

Fig. [5] shows the instability phase diagram of a single-layer 
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FIG. 5. (Color online) (a) The phase diagram of quasi-two- 
dimensional dipolar fermions in a single-layer configuration. The 
green dashed lines show the boundary between different multi- 
subband normal phases (No, Ni, . . .), the yellow shaded region is 
a density ordered phase and the thick lines on the N-DW boundary 
indicate the unstable wave- vector, q = 2A:i?,o (lower segment, sky 
blue), q = 2kF,i (middle segment, red) and q = 2A;f,2 (upper seg- 
ment, black), (b) The variation of unstable wave- vector along the 
N— DW boundary (black line). The blue, red and black dashed lines 
show twice theFermi momentum of the zeroth, first and second sub- 
band on the boundary. 



system as a function of rd and y/na^ . As we speculated be- 
fore, we find that the density-wave instability preceeds the rip- 
plon instability in the studied range of parameters. Therefore, 
the ripplon instability may only appear in the density ordered 
phase and form a zigzag phase (see Fig. [2]l). The plot next 
to the phase diagram in Fig. |5] shows the wave- vector of the 
unstable mode on the N-DW boundary. The switching of un- 
stable mode in the Ni can also be clearly seen: the density 
ordering wave- vector of a system in the Ni liquid phase is 
q = 2/ci?^o for \/na± < 0.25, however, it discontinuously 
jumps to q = 2kF,i for ^^a\_ > 0.25. The same behavior 
is observed in the N2 phase as well. We will investigate this 
behavior in Sec.lVl 



B. Instabilities of multi-layer systems 



As mentioned in Sec. II B the normal phase of multi-layer 
systems is very similar to single-layer systems in the IL limit, 
the only difference being existence of a mean-field shift of the 
subband energies due to direct inter-layer interactions. The 
inter-layer interactions, however, can drammatically affect the 
density wave fluctuations. In particular, one expects a more 
pronounced enhancement of both density wave and ripplon 
fluctuations. Analgous to the single-layer case, starting 
from the Nj phase, we again expect to see j + 1 peaks in 
the static density-density response functions at q = 2kF^o, 
q = 2/ci?,i, ' • Q = 2/ci?^j and a peak at q = kp^a + ^f,/3 
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FIG. 6. (Color online) The static density-density response function 
of a multi-layer system (^/nd = 1.25, Ni = 50) in the normal phase. 
Xdd = 2nh^Xdd/m and the blue and red planes denote q = 2A:f,o 
and q = 2A:f,o respectively, (a) No phase [vd = 1.265, y/na± = 
0.20] (b) Ni phase [r^ = 1.255, V^a± = 0.22] (c) Ni phase [r^ = 
0.845, y^a^ = 0.36]. In all plots, it is noticed that kz = modes 
experience the most softening. 

in the a/3-ripplon response function. The coupling between 
density- wave and ripplon modes at finite transverse momenta 
results in the mixing of these peaks such that traces of density 
wave peaks can be noticed in the ripplon response function 
and vice versa. In the following discussions, we keep the 
number of layers constant, Ni = 50, which is in the order of 
the typical number achievable in the experiments. 

Figs. [6] and [7] show the static density-density and 01 -ripplon 
response functions evaluated for three different points in the 
normal phase: the (a) plots correspond to a point in the Nq 
phase, the (b) plots are in Ni phase with a small population 
in the first excited subband and (c) plots are deep in the Ni 
phase. 

The plots in Fig. [6] indicate that the density- wave modes 
with zero transverse momenta experience most enhancement 
from the attractive inter-layer interactions. This is an expected 
result given that density-wave fluctuations are in-plane density 



FIG. 7. (Color online) The static 01 -ripplon response function of a 
multi-layer system (y^d = 1.25, Ni = 50) in the normal phase. 
Xoi = 27vh'^xoi/m and the green planes denote q = kF,o + kp,!- 
(a) No phase [rd = 1.265, y/na± = 0.20] (b) Ni phase [rd = 1.255, 
y^a± = 0.22] (c) Ni phase [rd = 0.845, y/na± = 0.36]. It is noticed 
that kzd = TV modes experience the maximum enhancements due 
to interactions. The peaks at n~^^^q ^ 3 for finite kz in plot (c) 
correspond to the softened density-waves which are coupled to the 
01-ripplons through parity-violating inter-layer interactions. 



modulations and at kz = 0, they are aligned across the lay- 
ers and thus, experience the maximum softening due to inter- 
layer attraction. We note that one expects the reverse scenario, 
i.e. maximum suppression of density-waves at kz = 0, had 
the inter-layer interactions been repulsive (such as multi-layer 
systems of 2DEG). 

On the other hand, the odd ripplons are expected to 
experience most softening at kzd = tt which corresponds to 
dimerization. At kzd = tt, the out-of-plane bumps of even 
numbered layers lie closest to the valleys of odd numbered 
layers, forming an energetically favorable configuration 
(shown schematically in Fig.|2|3). The slightly higher peak of 
01 -ripplon response function atkzd = tt compared to A:^ = 
is noticeable in Fig.|7})-c. The smaller peak in the 01 -ripplon 
response function (visible for < kz ^ it/ 2d) is due to 
coupling to the softened density- wave mode at q = 2kF,i- 
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FIG. 8. (Color online) The phase diagram of quasi-two-dimensional 
dipolar fermions in a multi-layer configuration (^/nd — 2, Ni = 
50). The black dashed line is the N-DW boundary in the single- 
layer configuration (refer to to caption of Fig.[5]for the description of 
the lines and symbols). 
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FIG. 9. (Color online) The phase diagram of quasi-two-dimensional 
dipolar fermions in a multi-layer configuration (y/nd = 1.5, Ni = 
50). The black dashed line is the N-DW boundary in the single-layer 
configuration (refer to to caption of Fig. |5] for the description of the 
lines and symbols). The hatched region is where the IL limit is not 
applicable (Ref. to Sec. [Till. 



In all of the studied cases, although the ripplon softening 
was found to be a more pronounced effect in multi-layer 
configurations compared to single-layer systems, the density- 
wave instability still precedes the ripplon instability. The first 
density- wave mode that becomes unstable has zero transverse 
momentum, implying that the density-wave and ripplon 
fluctuations are decoupled. Therefore, the density ordered 
phase to follow does not necessarily have out-of-plane order. 
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FIG. 10. (Color online) The phase diagram of quasi-two-dimensional 
dipolar fermions in a multi-layer configuration {sjnd — 1.25, Ni — 
50). The black dashed line is the N-DW boundary in the single-layer 
configuration (refer to to caption of Fig. [5] for the description of the 
lines and symbols). The hatched region is where the IL limit is not 
applicable (Ref. to Sec.|ll|. 



In the remainder of this section, we discuss the phase dia- 
grams of multi-layer systems for three inter-layer separations 
^/nd = 2, 1.5, 1.25. 

Fig. [8] shows the phase diagram of a multi-layer configu- 
ration with y/nd = 2 and Ni = 50. The dashed black 
line on the left plot indicates the N-DW boundary line of the 
single-layer system (copied from Fig. [5]). As mentioned be- 
fore, the first unstable mode is an in-plane density-wave mode 
with zero transverse momentum. We also find the most notice- 
able deviation of the N-DW phase boundary occurs for larger 
values of a±. The switching of the unstable wave- vector from 
q = 2/cf,o to g = 2/cf,i in the Ni phase is also found to oc- 
cur for a smaller value of ^Jna^ compared to the single-layer 
case, i.e. closer to the Nq-Ni boundary. 

Fig. [9] shows the phase diagram for ^Jnd = 1.5 and 
Ni = 50. The hatched region indicates the configurations 
at which the inter-layer tunneling is not negligible anymore 
and the approximation of independent layers is not justified. 
The physically interesting part of the phase diagram, however, 
lies outside of the hatched region. We notice that the N— DW 
boundary line deviates even further from that of single-layer 
systems. The switching point of the unstable wave- vector lies 
very close to Nq— Ni boundary. In other words, the Ni liquid 
phase always goes unstable due to the softened density-wave 
mode at g = 2A:f,i. Since =0 along the Nq— Ni transi- 
tion line, the unstable wave- vector can be arbitrarily small in 
the vicinity of the switching point (see Fig.[9]3). 

A more dramatic behavior is observed for smaller layer sep- 
shows the phase diagram for ^/nd = 1.25 
It is noticed that the N— DW boundary line 
becomes virtually tangent to the Nq— Ni transition line in the 
range 0.21 < ^/naj_ < 0.26. Along this part of the phase 
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Fig. 
= 50. 
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boundary, the transition to the inhomogeneous phase is driven 
by extremely long wavelength density- wave modes. 

In the next section, we approach the same problem again us- 
ing an approximate RPA-like formalism. Although we do not 
expect quantitatively reliable results, we still find that such 
an approach yields interesting analytical insights into some 
of the peculiar results of this section, in particular, the sud- 
den switching of the unstable mode along the N-DW bound- 
ary and the appearance of long wavelength unstable modes in 
multi-layer systems. 



V. INSIGHTS FROM THE RPA APPROXIMATION: 
NEGLECTING SHORT-RANGE EXCHANGE 
INTERACTIONS 

The major quantitative results of this paper was presented in 
the preceding section by numerically evaluating the response 
functions in the TDHF approximation. However, some of the 
results do not appeal to immediate intuition. In particular, 
(i) in single-layer systems, starting from the Ni phase, it is 
not clear why the the unstable density- wave abruptly switches 
from q = 2kF,o to q = 2/cf,i as the population of particles 
in the first subband is increased (see Fig. |5]) and, (ii) the ap- 
pearance of extremely long wavelength unstable density-wave 
modes along certain parts of the N-DW phase boundary in 
multi-layer systems is puzzling. In this section, we develop a 
simplistic and minimal model by applying successive approxi- 
mations to the TDHF formalism to derive an RPA-like expres- 
sion for the density-density response function, using which we 
will qualitatively explain the above findings. 

We start by noting that the main difficulty in obtaining ana- 
lytical expressions in the TDHF approximation is the exact 
treatment of exchange interactions, i.e. solving the Bethe- 
Salpeter integral equation. In the RPA approximation, on 
the other hand, one completely neglects the exchange in- 
teractions and this difficulty does not arise. However, the 
RPA approximation is not readily applicable to our problem, 
given that large cancellations are expected between the di- 
rect and exchange interactions of particle-hole pairs. This can 
be easily seen in the simplest case, i.e. a single-layer sys- 
tem in the two-dimensional limit (a± — > 0). In this limit, 
the only relevant interaction matrix element is Voo;Oo(q) = 
4V2^D^/3a^ - 27rD2^e-^l^l'^ + 0{D^q^a^). The ^-wave 
component of Voo;Oo(q) diverges in the limit a± 0. In 
a system of spinless fermions (which is the case here), the 
5'-wave interactions between the particles must vanish due to 
Fermi statistics and this cancellation only happens if one con- 
siders both direct and exchange interactions in a balanced way. 
This is clearly not the case in the RPA approximation. In these 
cases, it is customary to resort to heuristic methods to capture 
the exchange effects in an approximate way. Hubbard-type 
many-body local field approximations are widely used in the 
study of electron liquid | 28 | and have also been generalized 
to quasi-two-dimensional systems |29|. Such approximations, 
however, essentially aim at improving the long wavelength be- 
havior of the response functions. In our problem, we are in- 



terested in the response to density wave fluctuations at wave- 
lengths in the order of the inverse Fermi momentum. There- 
fore, the many-body local fields used for electronic systems 
are not readily applicable to our problem and must be modi- 
fied. 

Since we are only interested in qualitatively relevant re- 
sults in this section, we take the easiet route and argue that 
by simply removing the ^-wave component from all of the in- 
teraction matrix elements, the RPA formalism yields reason- 
ably decent values for the density-density response function 
3.tq^2 kpj- This claim can be justified by investigating the 
Bethe-Salpeter equation for the irreducible polarization with 
more care. For the clarity of argument, we consider the single- 
subband limit first, where the bookkeeping of subband indices 
can be obviated. Taking the static limit, Vn ^ 0, and defin- 
ing /(q,k) = n5o;oo(q.^O^;k)/n^o^.oo(q,iO+;k), one can 
rewrite Eq. ([37]) as: 

/(q,k) = l-| ^w(k'-k)n(0)(q,k')/(q,k'), (47) 

where ii(k) = Voo;Oo(k) and: 



n(0)(q,k)^n(°)oo(q,iO+;k) 



(ek.o) 



k+q,Oj 



^k,0 - ^'k+q,0 + ^0+ 



(48) 



For concreteness, we set q = According to Eq. ( [48] ) 

and as shown in Fis. is singular at ko = 

— kpx and we expect the most important contributions to the 



integral on the right hand side of Eq. ( [47] ) to result from the 
regions in the vicinity of ko. Thus, we may approximately re- 
place u{\i' — k) with u{\i' — ko) in the integrand. On the other 
hand, E^qjooI^^f:^, ^0+) = \, f(2kFX,\<i')Ii^^\2kFX,k') 
by definition, in which we may again approximately replace 
f{2kFX^'k') with /(2/ci?x,ko) according to same argument. 
Combining both approximations, we find that the final recipe 
is to replace ii(k' — k) with u{0) in Eq. (47 ), i.e. to keep only 



the long-range exchange interactions. We emphasize that the 
above argument is special to the analysis of |q| ^ 2kF modes. 

Once the short-range exchange interactions are neglected, 
the Bethe-Salpeter equation can be trivially solved. Combin- 



ing this results with the Dyson's equation, Eq. (40), we find 
that the only effect of the long-range exchange interactions is 
to remove the ^-wave component from the interaction matrix 
elements, as we expected. 

The above argument can be easily generalized to multi- 
subband and multi-layer systems using a matrix notation and 
we omit it here. In brief, we find that the general recipe is to 
simply make the substitution Vpzy;^o-(k — k') Vpiy;^cr(0) 
in the Bethe-Salpeter equation, yielding the following linear 
algebraic system of equations: 



(49) 



We have dropped the common arguments and subband indices 
for brevity in the above equation. Also, matrix multiplication 
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FIG. 11. Density plot of nQQ^QQ(2A;F,o^, k) showing the singular 
behavior at k = —kF,ox. 



is implied in each pair of subband indices. The approximate 
exchange interaction matrix, Vxch, is defined as: 



Combining Eqs. ([49]) and ([42|), we get: 



(50) 



n(fe.) = n(o) +n(o)^y(/c.) _y^^^^n^/c.)^ ^5^^ 

To ensure no violation of conservation laws, the short-range 
exchange interactions must also be neglected in the self- 
energy corrections. However, the long-range direct and 
exchange intra-layer exchange cancel each other. As men- 



tioned in Sec. II B the direct inter-layer interactions merely 
shift the zero-point subband energies by a small amount and 
for simplicity, one may neglect such corrections as well. 
Therefore, self-energy corrections can be neglected alto- 
gether. We refer to this approximation scheme as RPAns for 
brevity, with the ns subscript indicating the absence of ^-wave 
interaction terms. 

The important features of the phase diagrams presented in 
the previous section can be captured by keeping only the first 
two subbands. We also restrict the forthcoming analysis to 
kz = 0, given that such modes become unstable first. Under 
such assumptions, Eqs. (51 ) and ( [44] ) yield: 



Xdd{q) 



where: 



-1 



t(o) 



(0)/ 



detu{q) 



(52) 



detn(g) = 1 - Votoo(^)n^?(^) - V!i-M^)<\q) 
+ U^^,\q)U^^,\q) (Votoo(^)Vitii(^) - Votii(^)') • (53) 

The effective interaction matrix elements, V^^.^;^, are defined 
as: 



v: 



eff 



■ N,-l 

E 

.n=-Ni^l 



V^/3;c.a(0), (54) 



and the bare static intra-subband polarization, ni^2(q), can 
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FIG. 12. The effective interaction matrix element Voo;Oo(^) vs. q. 
The green (upper), blue (middle) and black (lower) lines correspond 
io Ni = 1, 3 and 200 respectively. In all cases, a±/d = 1/15. 
The inset plot shows the Ni = 200 case at small values of q for 
clarity. In the IL limit (a± <C d), one can classify the length scales 
into four categories according to the behavior of effective interac- 
tions, as indicated on the figure. Category (0): [q < L~^] length 
scales larger than L. The layered structure of the stack is invisible 
to density-wave fluctuations in this length scale. Since we have set 
kz = 0, the in-plane density waves are all aligned across the layers 
and collectively behave like a single density-wave with an effective 

(Ref. to Eq. 



58 1. In 



dipolar interaction strength of (2Ni 
other words, the whole stack behaves like a single two-dimensional 
layer; Category (I): [L~^ < q < d~^] length scales smaller than 
L and larger than inter-layer separation d. Density wave fluctua- 
tions in any given layer interacts with a fraction {qL)~^ of other 
layers, hence, resulting in a constant scale invariant effective inter- 
action -27tD'^{2Ni - l)q X (qL)-^ ^ -AivD'^/d; Category (II): 



[d- 



< 



q < a^^] length scales smaller than d and larger than a^. In 



this regime, the inter-layer interactions are exponentially small (see 
Eq. [56| and the density waves only interact within the layers. Cat- 
egory (III): [q < aj^] length scales smaller than a±. Each of the 
interaction matrix elements (Voo-oo^ Voo;ii^ etc.) assume different 
non-universal constant values in the order of /a±. 



be evaluated analytically in the absence of self-energy correc- 
tions: 



nL^2(q) 




(55) 

In the above equation, {k^^^} are the Fermi momenta of a 
non-interacting quasi-two-dimensional gas, as shown in Ta- 
bleU 

In this simplified approach, the single-layer and multi-layer 
systems are treated likewise. The multi-layer effects are in- 
cluded in the effective interactions. In other words, V^^.^;^ is 
the sum of all intra-layer and inter-layer interactions. Study- 
ing the behavior of the effective interactions is thus a key step 
in understanding the difference between the phase diagrams 
of single-layer and multi-layer systems. We focus on the be- 
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u(o) 



k 



(0) 



a± < 







7^ < a± < ^ 

/27rn — ^ V2TTn 



y^27rn — n/a\ 



of Fig. [12] justifying this identity for wavelengths longer than 
ax- The second line of Eq. ( [53] ) can be neglected in light of 
the this observation, yielding the following simple expression 
for detn(^): 



TABLE L The Fermi momenta of the first two subbands of a non- 
interacting quasi-two-dimensional gas. 



havior of Voo oo(^)' which is find to be qualitatively identical 
to the behavior of the rest of the involved interaction matrix 
elements, Vli.ii{q) and Vooai(^)- Expanding Eq. (A4) about 
= 0, we get: 



(56) 



detn(g)^l-Votoo(^rf(^) 



Vitii(^)nrAg), (59) 



(0). 



Intuitively, the above equation implies that the net density- 
wave enhancement is the algebraic sum of RPA-like density- 
wave enhancement of each subband. 



using which we find: 



N, = 1 



Votoo(9) 



-2TTD'q coth ( ^ j Ni = oo, ^^^^ 



3 



where the neglected terms are 0{q'^aj_). For future reference, 
it is also useful to study the behavior of Voo;oo(^) finite 
Ni , and for wavelengths longer than inter-layer separation d. 
In this limit, we find: 



Votoo(^) 



-27rD'^{2Ni 



L-'<q<d-K 



(58) 

Again, the neglected terms are 0{q'^a±). We remind that 
L = Nid is the transverse size of the stack. While the 
effective interaction in single-layer systems has a linear 
dependence on q in the regime q <C a^^, its behavior is 
very different for long wavelength modes in multi-layer 
shows Voo oo(^) ^ function of q for three 
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FIG. 13. Plotof detn(<?) vs. 
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-1/2 



q for three points in the normal 
phase (Ref. to Eq. [52] and [53] and the following text for details): (A) 
A^o phase [r^ = 0.40, y^a± = 0.30] (B) J\fi phase [r^ = 0.50, 
y/na± = 0.41] (C) Afi phase [r^ = 0.50, y/na± = 0.50]. Refer to 
Fig.[T4]to see the location of these three points in the phase diagram. 
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Fig. 
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systems. Fig. 

different num'Bef of layers, Ni = 1 (green), 3 (blue) and 
200 (black). In the IE limit (a± <C d), one can classify the 
length scales into four regimes according to the behavior 
of effective interactions. These regimes are indicated in 
Fig. ( [T2] ) and a brief description for each is provided in the 
caption. Consequently, one can categorize the density wave 
fluctuations according to the same length scale classification 
and as we will see shortly, this is a key step in interpreting the 
features of the obtained phase diagrams. 

We start the analysis of RPAns with the simpler case of 
single-layer systems. The stability of the normal phase can 
be determined by looking at the behavior of detn(^), which 
is the term appearing in the denominator of the RPAns ex- 
pression for the density-density response function. For small 
interaction strengths, detn(^) ^ 1. Upon increasing the in- 
teractions, detn(^) decreases and eventually, crosses zero at 
some q, signaling the appearance of a softened mode. 

Generally, we found that the approximate identity 
Votoo (^) Vit 11 {q) ^ Vot 11 {q? holds wdl for aU q. In par- 
ticular, all interaction matrix elements behave similarly in the 11^^^ (2kF^Q) ^ and at g = 2kF^Q, the density- wave enhance 



shows the plot of detn(^) for a point in Nq phase 
(A) and two points in Ni phase, with small and large popu- 
lation of the first subband (B and C respectively). In it no- 
ticed that in (A) and (B), the most softened mode (i.e., smaller 
detn) is q = ^kp^o, while in (C), q = 2kF^i is the most 
softened. The shift of the unstable mode from q = 2kF,o 
to q = 2/ci7^i in the Ni phase can be explained in light of 
Eqs. ( [59] ) and Eq. ( [55] ). In order to simply the discussion, we 
note that as long as g < 2/ci?,a, U.al{q) has a positive constant 
value and rapidly falls for q larger than 2k f^^- Thus, one only 
needs to monitor detn(^) for q = 2/ci?,o and q = 2kF,i, where 
the product of the effective interactions and the bare polariza- 
tions is be the largest. There are two possible scenarios in the 
Ni phase: 

Case I (kF,i <C /cf,o)' this case corresponds to dilute quasi- 
particles in the first excited subband and consequently, the 
effective interactions (which increase linearly with momen- 
tum) are weak at q = 2k f,!- Therefore, the sum of RPA-like 
enhancements resulting from both subbands at q = 2kF,i is 
smaller than the enhancement resulting mainly from the ze- 
roth subband at q = 2kFo (see Fig.[T3^). Since kFi <C kFo, 

(0)/ 



limit q 



< 



according to the remarks given in the caption ments are mainly due to the interactions in the zeroth subband. 
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Case II {kp^i ~ kp^o)'- this situation arises when there is a 
significant population in the first excited subband, i.e. deep in 
the Ni phase. The scenario is reversed in this case and the sum 
of enhancements resulting from both subbands at = 2/ci7^i is 
larger than the enhancement resulting mainly from the zeroth 
subband at g = 2A:f,o (see Fig.[T3p). 

It is not hard to see that the second scenario may only 
happen if the rise of interactions is slower than the fall of 
density of particle-hole excitations as a function of q. The 
linear momentum dependence of dipolar interactions and the 
rapid fall of Ii^^^{q) fox q > 2kF,i guarantees the realization 
of this situation for large enough values of kp,!- 



Fig. [14] shows the approximate phase diagram of a single- 
layer system calculated in the RPAns approximation. The 
flatness of Nq-Ni and N1-N2 boundaries is due to ignoring 
the self-energy corrections in the normal phase. There is a 
striking similarity between this phase diagram and the one 
obtained by exact numerical calculation of TDHF response 
functions (Fig. [5|. However, the predicted value for the DW 
instability at a_\_ 0, r^^^ ^ 0.15, is more than a factor 
of two smaller than the same value predicted within TDHF, 
r™^ ^ 0.39. 



Fig. [T5] shows the approximate phase diagram of three 
multi-layer systems with different inter-layer separations ob- 
tained using the RPAns scheme. It is noticed that the non- 
trivial features of multi-layer phase diagrams, i.e. (1) indiffer- 
ence of N-DW boundary line to existence of multiple layers 
deep in the Nq phase, and (2) enhancement of density wave 
instability along parts of Nq-Ni phase for smaller inter-layer 
separations, are also present in the picture that RPAns sug- 
gests. 

The former feature can be explained by first noting that the 
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FIG. 14. The approximate phase diagram of quasi-two-dimensional 
dipolar fermions in a single-layer configuration in the RPAns ap- 
proximation. The dashed lines show the smallest rd for which the 
density- wave mode at ^ = 2kF,o or q = 2/cf,i become unstable. 
The pink line indicate the first unstable mode. The switching of un- 
stable density- wave mode in the Ni phase is noticeable. Refer to 
Fig. 13 for a plot of detn(^) for the three points marked in the dia- 
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FIG. 15. The approximate phase diagram of quasi-two-dimensional 
dipolar fermions in multi-layer configurations {Ni — 20) in the 
RPAns approximation: (a) yjnd — 3 (b) yjnd — 2 (c) yjnd =1.5. 
The dashed lines show the smallest rd for which the density-wave 
mode at ^ = 2A;f,o or ^ = 2A;f,i become unstable. The pink lines 
indicate the first unstable mode. 



studied range of inter-layer separations is such that ^/nd ~ 
0(1). Therefore, the unstable vector in the Nq phase, q = 
2A^F,o ^ 2V47rn, is almost an order of magnitude larger than 
d~^ and belong to category (II). The inter-layer interactions 
are irrelevant in this regime and the physics is identical to that 
of a single-layer system. 

As a side note, we remark that in order to see the effect of 
inter-layer interaction in the single- subband limit (Nq phase), 
one must choose d such that ^/nd <C 1. In particular, in the 
limit ^/nd <C A/"^"^, the unstable modes will lie in category 

(0) and the density wave instability will be driven by fluctua- 
tions whose length scale is larger than the transverse length of 
the stack. As mentioned earlier, the effective interactions are 
enhanced proportionally to the number of layers in this limit 
and as a result, we expect the interaction strength required for 
the onset of density- wave instability to be reduced by a factor 
of^N'K 

The latter feature, i.e. appearance of long wavelength un- 
stable modes close to Nq-Ni boundary can be explained as 
follows. We discuss the simpler case of A'^ ^00 first, in 
which all q ^ lie inside category (I), i.e. where the effec- 
tive interactions assume a constant value of —AirD'^/d. Ex- 
istence of a small particle density ni in the first excited sub- 
band will result in the appearance of long wavelength gapless 
particle-hole excitations. The length scale associated to these 
modes can be very large and may as well lie within category 

(1) for small enough ni, i.e. q = 2A:fi ~ 2v'47rni ^ d~^. 
Since the density of long wavelength excitations is finite in 
two dimensions, i.e. lim|q|^o n'^^^(q) ~ 0{m/2iTh^), they 

(0)/ 



gram. 



will have a finite RPA-like contribution of IIii (2fe/ ?.T ) 
Vff.ii(2/cF,i) ^ 2mD^/dh^ to detn (see Eq. 



59). For 



small inter-layer separations, this contribution can oe large 
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and result in density-wave instability. 

In the limit Ni ^ oo, these modes appear exactly along the 
Nq-Ni boundary where q = 2kF,i = 0. The largest layer sep- 
aration, (imax, for which such long wavelength modes appear 
can be easily determined. At d = d^g,^, the = unstable 
mode appears only at one point, viz. at the intersection of Nq- 
Ni and N-DW lines. Therefore, both q = 2kF,o and q = 
are unstable at this point (see Fig. 15 i). The RPAns instability 
condition at g = yields: 



detn(0+) = 1 



= 1 



0, 



^max 



(n^?(0) + nff(0)) 



m 



^max Tt/i^ 



(60) 



and the instability of q 

detn(2/cF,o) 



= i-v; 

= 0. 



00:00 



2/cf,o = 2v47m yields: 



27rn2 



(61) 



In the above equation, the effective interaction must be eval- 
uated on the Nq-Ni boundary, i.e. aj_ = 1/V27rn, i.e. The 
simultaneous solution of these equations yields: 



Vndn 



2.209, Td - 0.5523. 



(62) 



For d < (imax, the = unstable modes appear along a finite 
interval on the Nq-Ni boundary (see Fig. 10 and Fig. 15 3-c). 

The prediction of (imax within RPAns is significantly larger 
than the one inferred from TDHF calculations presented 
earlier (see Fig. [9j ^/nd™^^ 1.5 and rd o:^ 1.35). This 
deviation is again due to approximate treatment of exchange 
interactions. 

We conclude this section by briefly studying the scaling 
dependence of the wavevector of the long wavelength unsta- 
ble modes discussed above on Ni. For finite Ni, q = 2/ci?,i 
mode lies inside category (0) if the first excited subband is 
dilute enough. In this limit, the whole stack behaves collec- 
tively like a single pancake, with an effective interaction of 

-27r(2A^/ - l)D'^q. Assuming q < the RPAns insta- 
bility condition yields: 

detn(g) ^ 1 - 2 X ^ 27rD\2Ni - l)q = 0. (63) 

Solving for q, we get: 



^/nd 



2{2Ni-l)rd' 



<rd< 0.55. 



(64) 



The constraints imposed on rd in the above equation result 
from two requirements: on one hand, the solution must satisfy 
< On the other hand, the required for instability of 
this mode must be smaller than that required for the instability 
of the q = 2/ci? Q mode, which is ~ 0.55 in the vicinity of 
the Nq-Ni boundary and for d not much less than (imax (see 



Eq. 62 and Fig. 15 ). In the limit Ni ^ oo, the unstable wave- 
vector becomes 0. 



VI. DISCUSSION AND CONCLUSIONS 

In this paper, we studied the mean-field density order- 
ing instabilities of quasi-two-dimensional fermionic polar 
molecules in both single-layer and multi-layer configurations. 
The dipole moments of the molecules were assumed to be 
aligned perpendicular to the confining planes using a dc elec- 
tric field. We located the instabilities by evaluating various 
linear response functions in the liquid phase and searching 
for the softened modes. We considered both in-plane and 
out-of-plane density ordering instabilities, as schematically 
depicted in Fig. [2] 

In all of the studied cases, the instability of the in-plane 
density wave modes was found to precede the instability of 
out-of-plane "ripplon" modes, although the latter modes were 
also softened to some degree. We also found that leading 
unstable mode in multi-layer systems has a zero transverse 
momentum, i.e. the in-plane density-waves are aligned across 
the layers. 

In multi-layer configurations, an interesting finding was the 
enhanced density wave instability driven by dilute quasipar- 
ticles of the first excited subband. By analyzing the effective 
interactions at various length scales in Sec. |Vj we found that 
these dynamical instabilities are associated to the softening 
of low-energy particle-hole excitations whose wavelength 
is comparable to or larger than the transverse size of the 
system, L. On one hand, the density of such excitations is 
finite due to the underlying two-dimensionality of the system. 
On the other hand, their effective interaction is enhanced 
proportionally to the number of layers due to their long 
wavelength. Hence, they produce a significant density-wave 
softening effect. 

Another interesting feature of the phase diagram of 
both single-layer and multi-layer configurations is the non- 
monotonicity of the N-DW phase boundary as a function of 
transverse confinement width (see Figs. [5j [8pQ|). A phase 
diagram with similar qualitative features had been predicted 
before for quasi-2DEG using density functional theory with 
Perdew-Zunger type exchange-correlation energy IJTII . Thus, 
we expect this feature of the presented instability diagram to 
persist in the true phase diagram, i.e. when the correlation ef- 
fects are also taken into account. The protrusion of the liquid 
phase inside the density ordered phases allows realization of 
the following interesting experimental scenario: starting from 
a point such as rd = 1.45, y/na^ = 0.15 in a density ordered 
phase of a single-layer system, the homogeneous liquid phase 
can be revived upon relaxing the trap. However, the liquid 
state becomes unstable again upon relaxing the trap further 
(i.e. by crossing the red boundary line in Fig. |5]) toward a 
different density ordered phase. 

As mentioned in the introduction, Yamaguchi et ah 1 14 1 and 
Sun et ah 1 15 1 have recently studied the density- wave instabil- 
ity of a (strictly) two-dimensional system of polar molecules 
(a^ 0) within the RPA approximation. Their study, how- 
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ever, includes the more general case of tilted dipoles and fi- 
nite temperatures. At zero temperature and dipole tilt angle, 
their results indicated that the density-wave instability occurs 
at Td ~ 0.17. The RPA-like approximation used in Sec. [V| 
yields ^ 0.15 in the limit a± ^ which is in good agree- 
ment with the result of Ref. 1 14 |. 



The TDHF results presented in Sec. IV however, predicts 
Td ~ 0.39 in the same limit which is more than a factor of 
two larger than the prediction of the RR\ approximation. 
Hence, we conclude that an exact treatment of short-range 
exchange interactions is important for quantitatively reliable 
predictions of phase transitions in dipolar systems. 

The instabilities predicted in the mean-field picture must be 
interpreted with care. On one hand, one must consider the 
possibility that the actual phase transition is first-order. In this 
case, the mode softening criterion does not indicate the true 
transition but signifies the spinodal line (i.e. the end of liquid 
metastability region) and the actual phase transition occurs for 
smaller values of rd- Typically, the spinodal line lies close to 
the actual transition line |38 | and therefore, we do not expect 
the above issue to be a major source of error in the obtained 
phase diagrams. 

The main shortcoming of the present analysis lies in the 
mean-field approximation and absence of correlation effects 
in the liquid phase. It is known that the mean-field mode 
softening analysis often predicts the transition to the symme- 
try broken phases too early. For instance, the Wigner crystal 
phase of 2D electrons with 1/r Coulomb interactions is found 
to become stable for > 1.44 in the mean-field approxima- 
tion 1301 while a more realistic quantum Monte-Carlo calcu- 
lation yields a value of > 38 |31|. Thus, we expect that 
the instability lines shown in Figs. [5] and 8pQ will be shifted 
to larger values of rd once correlation effects are taken into 
account. Since mean-field predictions improve by increasing 
the dimensionality, this correction is expected to be smaller in 
multi-layer systems compared to single-layer systems. Never- 
theless, we expect that the mean-field transition lines obtained 
here will describe sharp crossovers to the regime of strong 
short range crystal correlations (with no long-range order) in 
real systems, with the actual phase transitions following for 
larger values of rd- 

We remark that the regime of strong short-range crystal cor- 
relations with no long-range order is reminiscent of the pseu- 
dogap phase of fermions with strong attractive interactions. In 
the latter case, one finds short-range pairing correlations but 
no true long-range order, i.e. no condensation of molecules. 

While the full analysis of such strongly correlated "pre- 
formed density- wave" state is outside the scope of this paper, 
the mean-field analysis presented here is an indispensable 
first step toward the study of this intriguing state. 

At the time this paper is being written, the experimental 
verification of the presented results can still be challenging. 
The maximum dipolar interaction strength accessible in the 
experiments is rd ~ 0.05, which belongs the experiments of 
the group at JILA with KRb polar molecules. Observation 
of density-wave instability either requires production of 



denser gases or using molecules with larger dipole moments 
(the permanent dipole moment of KRb is 0.55 D). Further 
progress in experiments with LiCs |0 and RbCs (391 polar 
molecules whose permanent dipole moments are 5.5 D 
and 1.25 D respectively, are among the other promising 
experimental directions toward observation of the effects 
discussed in this paper. 



We remark that in the same way ultracold atoms were 
utilized as a simulator for confined electrons with effective 
short-range interactions and shed light on the Hubbard 
model, ultracold polar molecules may be utilized as a tool to 
address unsettled questions regarding the nature of transitions 
to density ordered phases, intermediate strongly correlated 
states (such as the electron nematic phase) and microemulsion 
phases (such as stripes and bubbles) 1341 . 



Throughout this work, we had assumed that the stable 
phase in the weakly interacting regime is the normal liq- 
uid phase. Neglecting the weak high angular momentum 
superconducting phases which may only appear at very 
low temperatures, this assumption is valid for single-layer 
systems. In multi-layer configurations, however, the normal 
phase is known to be unstable toward dimerized pseudo-gap 
and inter-layer superfluid phases ll25ti27l at low temperatures. 
The passage through these phases occur through Ising-like, 
and Berezinskii-Kosterlitz-Thouless phase transitions, re- 
spectively. In our study, the configuration which is most likely 
to be in a superfluid phase in all of the presented multi-layer 
phase diagrams (Figs. [SpO} excluding the hatched regions) 
is ^/nd = 1.25, rd ~ 0.8 and y/na± ^ 0.28. For such a 
configuration, the critical temperature for BCS transition is 
estimated to be Tc/T^ ^ 0.013 |40| using the results of 
Ref. 1251 . Therefore, the temperature chosen in this study, 
T/Tp^ = 0.02, is above the inter-layer pairing transition 
and our assumption about the stability of the liquid phase for 
weak interactions is justified. 



The competition between inter-layer pairing instability and 
density-ordering instabilities at lower temperatures, or for 
systems with smaller inter-layer separations, is an interesting 
topic for future works. The results presented here can also be 
easily extended to tilted dipoles. Reducing the intra- subband 
repulsion and enhancing the inter-subband repulsion, tilting 
the dipoles may reverse the order of density-wave and 
ripplon instabilities. The competition between intra-layer 
/?-wave superfluidity which is expected to occur for large tilt 
angles |41 1, ripplon and density- wave instabilities is another 
interesting topic for future research. 



Note added: After the completion of this work, we became 
aware of a recent related work by Zinner et al | 42 1 in which 
they study the density-wave instability of stacks of strictly 
two-dimensional polar molecules (a± 0) within the RPA 
approximation. 
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Eq. ( |A4| ) is a simple Gaussian integral while the second dou- 
ble integration can be easily evaluated by changing variables 
io T] = {z ^ z') /2 and ^ = {z — z')/2 and successive integra- 
tion by parts. We just quote the final result here: 

l/oo;Ooiq, — — e ^ 



3a I 



■^I)2|q|F+(|q|ax,V«±), (A5) 



Appendix A: Analytical expressions for the interaction matrix 
elements 



In this appendix, we provide analytical expressions for 
inter-layer interaction between quasiparticles in the first two 
subbands. The interaction between quasiparticles in higher 
subbands may also be calculated using the same method. 

Using Eq. ^ and approximating the Wannier's wavefunc- 
tions by shifted harmonic oscillator wavefunctions, the effec- 
tive interaction between particles confined to planes z = 
and z = I can be expressed as: 

Vc./3;7A(q;0 = if ^zMz)(l)(3{z)e-'^^ 



J dz' (l)^{z' - I) (})x{z' - I) i 



Akz 



(Al) 



where (j)a{z) is a'th harmonic oscillator wavefunction and: 



'Kiip(^,q) = j dz& 



(A2) 



Note that V^'^.'^^(q) = V^^.^x (q; (n - m)d). Evaluating the 
k integral, we get: 



(A3) 



Plugging Eq. (A3 ) into Eq. ( Al ), we get: 



Va/3;7A(q;0 = 

j dz (j)oc{z)(j)p{z)(j)^{z -l)(j)x{z -I) 

-2nD^\ci\ j dz j dz' e-\'^\\'-''\(^^{z)(^^{z) 

X (l)^{z' -l)(l)x{z' -I). (A4) 

At this point, one may proceed by finding a generating 
function for Va/3;7A(q;0 through expressing the Hermite's 
functions appearing in the harmonic oscillator wavefunctions 
in terms of their generating functions 1 18 1. Since we are inter- 
ested in the first few matrix elements here, we find it is easier 
to evaluate the required integrals directly. The first integral in 



+ y27rL>2a_L|q|2F_(|q|a_L,Va-L), (A6) 



Voo;oi(q;0 = — '-^ — e 



Voo;ii(q;0 = ^^^^e-'W(i + ;/„^) 
+ ^|q||2v^|q|a^e-'V2ai 

-7r(2+|q|2ai)F+(|q|ai,Va±), (A7) 



Voi;ii(q;0 



3al 
0FZ)2|q|2 



q| < — 4l/aj_ e 



Vii;ii(q;0 



4a_L 

M2+|q|'ai)J^-(|q|a±,Va±) \, (A8) 



q| <! - 4|q|aie-' /^^l (3 + |q|^a^ + t^/a^ 



4 V 2 



27r(2+|q|X)^+(|q|a±,Va±) )■ (A9) 



In the above equations, F±{x, y) is defined as: 



F±{x,y) = e 



V/2 



e(-.)V2Erfcf^ 
\/2 



±e 



(-+.)V2Erfc('^ 
V2 



. (AlO) 



Appendix B: Numerical Solution of the Bethe-Salpeter Equation 

A major difficulty in evaluating response functions in the 
TDHF approximation is solving the Bethe-Salpeter equation 
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FIG. 16. Adaptively generated grids for integrations involving (a) 
^oo;7A, (b) hoi-^x, (c) hiQ-^x and (d) hn-^x for a system in the Ni 
phase [rd = 1.35, y^a± — 0.25] and ^ 4.78. ^cc(y) = 



resulting from the ladder diagram summations, Eq. ( 37 ), in or- 
der to obtain the irreducible polarization H*^.^;^. The book- 
keeping of subband indices in quasi-two-dimensional systems 
is an additional difficulty. Nevertheless, the problem is essen- 
tially a system of coupled Fredholm integral equations of the 
second kind, which can be efficiently solved using numerical 
methods such as the Nystrom method f43l. In this method, 
one approximates the integrations with using quadrature for- 
mulas and solves the resulting (large) system of linear equa- 
tions. 

In the approach used in this study for locating softened 
modes, one is only interested in the response functions in the 
static limit, i.e. ^ 0+. In this limit, the bare polarizations 



appearing in Eq. ( 37 ) will have integrable singularities at the 
intersection of particle and hole Fermi surfaces. For exam- 
shows the single singularity at ko = —kp^ox for 



11 



pie. Fig. 

q = 2kF,ox. The dressed polarizations may have additional 
singularities associated to the softened collective modes. The 
single-particle and collective singularities can be separated by 
simply dividing the irreducible polarization by bare polariza- 
tion. We define: 



^a/3;7A(q,^^n;k) = 

(n(0))"' (q,z^n;k)n;,,.^;,(q,z^,;k), (Bl) 



where the inverse of the bare polarization is de- 
fined only with respect to subband indices, i.e. 



Recasting Eq. (37 ) in terms of hafd^^x, we get: 



^l) = ^a-f^^X 



I 



V,ft«^(k'-ki) 



,(q, iVn] k') /lpa;7A(q, ^^n] k')- (B2) 

In the absence of interactions, hafd^^x is simply the identity 
operator in the space of subband indices. Since the bare 
polarization appears in the integrand, the integrable singular- 
ities associated to gapless particle-hole excitations will not 
result in any singularity in haf3;^x- On the other hand, if the 
Fredholm determinant of the above integral equation vanishes 
at some q, i.e. det [l + Vn*^^^] = 0, h^^-^x will be singular 
at that q. In fact, this condition can be used as a practical 
criterion for locating the softened collective modes. Thus, the 
single-particle poles are absent in ha^-^yx and it effectively 
represents the many-body correctiosn to the bare polarization. 

For non-smooth integral kernels, as it is the case here, rapid 
convergence of Nystrom mehod is only achieved if one em- 
ploys adaptively generated integration quadratures that prop- 
erly handle the integrable singularities and fast variations of 
the integral kernels. The singular points must be avoided and 
a finer mesh must be used in the proximity of the singulari- 
ties and sharp variations of the integrand. We implemented 
the adaptive mesh refinement (AMR) algorithm described in 
Ref. [? ] on a square-based mesh to generate the integration 
quadrature. For each q, a uniform rectangular grid was gener- 
ated and adaptively refined until the relative integration error 
was smaller than 10~^. One may generate a single "global" 
quadrature that handles the irregularities of the various inte- 
gral kernels appearing in Eq. ( |B2| ) corresponding to different 
choices of subband indices. However, a more efficient ap- 
proach can be devised by utilizing the parity conserving nature 
of intra-layer interactions. For instance, when only the first 
two subbands are relevant, there is no subband hybridization 

and ^alnx ^ 



^ ^a-i^^x^^al- Therefore, hp^r-^x only appears 
in conjuction with Hp^ in the intergand of Eq. (B2) and con- 



sequently, one may produce four separate integration quadra- 



tures for h[ 



00;7A. ^01;7A. ^10;7A 



and h 



11;7A; 



each of which 



has about half the number of points of a globally applicable 
quadrature. Fig. [16] shows an instance of the adaptive grid 
generated in this fashion. 

In all of the studied cases, the algorithm produced a mesh 
containing ~ 5000 (or less) points before the stopping crite- 



ria was fulfilled. The integrals appearing in Eq. (B2 ) was then 
approximated using the generated quadrature and reduced to a 
linear system. The linear system was solved using LU decom- 
position. Once hais^jx was calculated, the irreducible polar- 
ization diagrams were finally evaluated by multiplying h^fs^^x 
by the bare polarization and summing over ki : 



/d k 
^^n2];p^(q,^^n;ki) 

X /ip,.;7A(q,^^n;ki). (B3) 



The previous generated quadratures can be utilized to evaluate 
the above integral as well. 
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